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ABSTRACT 

NASA Glenn Research Center (GRC) is funding Cleveland 
State University (CSU) to develop a reliable Computational 
Fluid Dynamics (CFD) code that can predict engine 
performance with the goal of significant improvements in 
accuracy when compared to one-dimensional (1-D) design 
code predictions. The funding also includes conducting code 
validation experiments at both the University of Minnesota 
(UMN) and CSU. In this paper a brief description of the work-in- 
progress is provided in the two areas (CFD & Experiments). 
Also, previous test results are compared with computational 
data obtained using: 1) a 2-D CFD code obtained from Dr. 
Georg Scheuerer and further developed at CSU and 2) a multi- 
dimensional commercial code CFD-ACE+. The test data and 
computational results are for: 1) a gas spring and 2) a single 
piston/cylinder with attached annular heat exchanger. The 
comparisons among the codes are discussed. The paper also 
discusses plans for conducting code validation experiments at 
CSU & UMN. 

INTRODUCTION 

The type of flows that are normally encountered in Stirling 
Engines are: 1) Unsteady (Oscillatory Flow/ Oscillatory 

Pressure), 2) 2-D/3-D Flows (Complex Geometry), 3) 

Compressible (low Mach Number), Laminar/Transition/ 
Turbulent, 4) Unsteady Conduction/Convection Heat Transfer, 
5) Sudden changes in cross-section, 6) Isothermal/adiabatic 
boundary conditions, and 7) Single phase with no chemical 
reaction. On the other hand the losses encountered in Stirling 
engines (Tew 1988 & Tew and Geng 1992) are: 1) Gas spring 
and working space hysteresis losses, 2) Gas-to-heat exchanger 
heat transfer inefficiencies, 3) Viscous (pressure drop) losses, 
4) Appendix gap losses (shuttle and pumping), 5) Mixing losses 
due to mixing of gases at different temperatures, 6) Conduction 
losses from the hot- to the cold-end of the engine. In addition, 
there are losses due to differences in flow distribution: a) from 
one heat exchanger flow passage to the other and/or b) across 
the regenerator (1-D codes assume uniform flow). 

A successful CFD code that simulates the engine 
adequately should capture the characteristics of the flow and 
the type of losses described above. 


Stirling One-Dimensional Codes 

During the late 70’s and early 80’s, NASA developed it’s 
own finite-difference Stirling engine performance code for use in 
monitoring the work of its Stirling engine contractors (Tew, et. 
al.1978 & Tew 1983). Later NASA gained unlimited rights to the 
MTI developed HFAST (Huang 1993) harmonic code, via a 
contract with MTI. NASA also purchased the GUMPS Stirling 
engine code from David Gedeon. GLIMPS (now developed into 
Sage, Gedeon 1995) has been the primary design tool used by 
the Stirling Technology Co. (STC) in recent years. HFAST and 
GUMPS (now Sage) were both more time-efficient and user- 
friendly than the NASA’s finite-difference code; and they were 
both being used to design real hardware. 

Tew and Geng, 1992 showed comparisons of the losses 
calculated by GLIMPS and HFAST for the 12.5 kWe 
Component Test Power Converter (CTPC). Although the overall 
engine power and efficiency predictions were quite close, there 
were substantial differences in some of the loss calculations. 
For example GLIMPS calculations suggested that cylinder 
hysteresis power losses were about 10% of the indicated power 
while HFAST calculated a much smaller value (~3%) for this 
loss. HFAST predicted larger viscous and mixing losses than 
GLIMPS, and so the overall engine performance predictions 
were similar. Geng and Tew, 1992 showed additional 
HFAST/GLIMPS comparisons, both “calibrated" and 
“uncalibrated”, for the 1.2 kW indicated power RE-1000 engine 
and the 14 kW indicated power Space Power Demonstrator 
Engine (SPDE). GLIMPS also calculated that cylinder 
hysteresis losses were about 10% of the indicated power for 
these engines. 

The Sage commercial code (Gedeon 1995) was introduced 
about six years ago. The Sage Stirling-cycle modeling software 
is the latest in a line of commercial software developed by 
Gedeon Associates. It is a direct descendant of the GLIMPS 
software, which was used widely within the Stirling industry for 
nearly ten years. Sage introduced a drag-and-drop visual 
interface where a user could assemble complete machines from 
standard components, such as pistons, cylinders, heat- 
exchangers, etc. Sage also introduced an interactive 
optimization capability built into the visual interface. 
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CAST and Modified CAST Code 

A CFD code, CAST, received by Dr. Mounir Ibrahim from 
Dr. Georg Scheuerer, was developed by Peric and Scheuerer 
1989. CAST, an acronym for Computer Aided Simulation of 
Turbulent Flows, is a computer program which uses the finite 
volume method for predicting two-dimensional flow and heat 
transfer phenomena. The CAST code is written in FORTRAN 
IV. CAST is similar in structure to existing fluid flow prediction 
procedures like TEACH (Gosman and Ideriah, 1976) and TEAM 
(Huang and Lewchziner, 1983); it differs from those codes in 
the following ways (1) it has a co-located variable arrangement, 
(2) it has a different discretization scheme, (3) the solution 
algorithms are different for the linear equation systems arising 
from the discretization, and (4) the pressure-velocity coupling is 
adapted to the co-located variable storage. CAST solves the 
Navier-Stokes equations. For turbulent flows, the Reynolds- 
averaged Navier-Stokes equations are solved in connection 
with the k-e, two-equation turbulence model of Launder and 
Spalding. 

The CAST code was used by several Cleveland State 
University (CSU) graduate students working under the direction 
of Dr. Mounir Ibrahim for two-dimensional modeling of “Stirling 
machine like" components. This work was conducted under 
grants from NASA Glenn Research Center (which was NASA 
Lewis at the time of the work) in the period from 1989 to 1994. 

Below is a summary of the code development topics: 

1) Oscillating inlet & outlet velocity conditions, 2) Sudden 
change in channel cross section, i.e. sudden expansion in one 
half of the cycle and sudden contraction in the other half, 3) 
Different Low-Reynolds Number, k-e, turbulent models, 4) A 
model for laminar-transition-turbuient in a pipe flow, 5) 
Incompressible-flow with moving boundaries and, finally, 6) 
Compressible-flow with moving boundaries. 

This work is summarized below: 

Ahn 1992, added a low-Reynolds-Number k-e turbulence model 
to CAST for purposes of conducting studies of laminar to 
turbulent transition in pipe flows. This modified CAST code was 
used in cooperation with experimental studies at the University 
of Minnesota to help determine that upon relaminarization of 
turbulent pipe flow over a periodic cycle, the laminar flow at the 
beginning of the new cycle had attained the characteristics of a 
steady, uniform velocity distribution (in the radial direction). 

Hashim, 1992, configured the CAST code to simulate 
oscillating flow and heat transfer in parallel-plate channels with 
a sudden change in cross-section. The flow was assumed to 
be laminar and incompressible with the inflow velocity uniform 
over the channel cross-section but varying sinusoidally with 
time. Ibrahim et. al., 1992, using the computational results, 
were able to determine instantaneous friction factors and heat 
transfer coefficients for laminar oscillating flow between parallel 
plates with a sudden change in cross-section. It was found that 
instantaneous friction factors and heat transfer coefficients 
deviated substantially from the steady-state values for the same 
dimensionless flow parameters. 

Kannapareddy, 1993, used CAST to model laminar, 
incompressible, oscillatory flow in the heater, regenerator and 
cooler of the NASA Stirling Space Power Research Engine 
(SPRE). The heater and cooler tubes were modeled as circular 
pipes with isothermal walls; the regenerator was modeled as 
two-parallel plates. Although the flow was oscillatory, when the 
flow was into the heater or cooler, the input temperatures were 
assumed to be constant. The study was carried out for a wide 
range of Maximum Reynolds numbers (based on maximum flow 


velocity), Valensi numbers (non-dimensional frequencies), and 
relative amplitudes of fluid displacement in the component of 
interest; the ranges chosen were based on the operating 
ranges of the NASA SPRE. The instantaneous friction factor, 
wall heat flux and heat transfer coefficient were examined. It 
was concluded that the friction factor and heat transfer 
coefficients are larger under oscillatory flow conditions for larger 
Valensi numbers (i.e., larger non-dimensional frequencies). 
Also the thermal efficiency of the heat exchangers decreased 
for the lower fluid displacement values and the steady-state 
definition for the heat transfer coefficient did not appear valid for 
use with oscillatory-flow ranges studied. A concise presentation 
of this work was reported by Ibrahim and Kannapareddy, 1992. 
In work that was not documented, Kannapareddy modified 
CAST to incorporate a moving boundary model to allow 
simulation of piston motion. These modifications were based 
upon information in Ferzigerand Peric, 1997. 

Bauer, 1993, used CAST with the low-Reynolds number 
k-e, turbulence model option (Lam-Bremhorst form), to assist 
Ibrahim et al., 1994, in developing an empirical transition model 
that could be used to predict when in the University of 
Minnesota (UMN) periodic-pipe-flow experiments (Qiu and 
Simon, 1994) transition laminar to turbulent flow would occur. 
The empirical transition model was used to activate the 
turbulence model at the appropriate time within the cycle for a 
given axial location in the tube. This analytical technique for 
modeling unsteady flow and heat transfer in Stirling engine 
heater and cooler tubes (Simon, et al. 1992, and Ibrahim, et al., 
1994) is now an essential component of the 1-D system 
simulation code Sage (at that time called GLIMPS). The tests 
used for this development showed important features of 
oscillatory flow, which led to enhanced understanding of the 
mechanisms for losses in Stirling systems. Such features 
include transition from the laminar flow to turbulence in each 
cycle of oscillation and the rapid dissipation of the turbulence, 
which remains from one cycle as a result of the extremely 
strong axial straining of the fluid upon temporal acceleration at 
the beginning of the following cycle. This collaborative effort 
showed the value of combining computational results with 
experimental finding for better understanding of Stirling engine 
flow and heat transfer processes. 

In searching for a compressible flow simulation techniques 
suitable for simulation of Stirling machines, the literature was 
first searched for papers relating to low Mach number 
compressible flow.” Many such papers were obtained. For 
example, papers by Pletcher and Chen, 1993, Chen and 
Pletcher, 1991, Sesterhenn, Muller et al.1993, Turkel, 1987, 
Weiss and Smith, 1993, and Horibata, 1992 were reviewed! 
These papers appeared to take the approach of starting with a 
fully compressible model (including acoustics) and making 
suitable modifications to the solution technique so that the 
calculations could be extended to low Mach number flow. If any 
of these approaches had been taken to modify the CAST code, 
it appeared that it would be almost equivalent to starting over 
and developing a new code. And, these approaches did not 
appear to offer the time benefits that could be achieved by 
complete elimination of acoustic phenomena. 

For ail practical Stirling engines with which these authors 
have familiarity, flow occurs at very low Mach numbers (usually 
much less than 0.1). However, compressibility is very important 
because volume changes due to piston motion in an enclosed 
volume produce large changes in the engine pressure level and 
density over the cycle. Due to the compact nature of Stirling 
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machines and the moderate frequencies of operation 
(< 125 Hz), spatial pressure variations at any given time are 
relatively small so that spatial density variations are due almost 
entirely to spatial temperature variations. Calculations have 
shown that it appears to be a good approximation to assume 
that the speed of sound is infinitely fast inside most Stirling 
machines. Some authors, such as Organ, 1992, appear to 
contest this assumption 

The literature search continued for a compressible flow 
approach based on the assumption that acoustics are 
negligible. An excellent paper explaining such an approach 
was found in the 1997 Journal of Fluid Mechanics entitled “A 
model of unsteady subsonic flow with acoustics excluded" by 
AT. Fedorchenko, 1997. The derivations include a set of 
equations for simulating subsonic flow of a heat-conducting, 
viscous gas. The simplifications in the Navier-Stokes equations 
used to eliminate acoustics include: (1) Pressure at any time 
and spatial location is split into a mean pressure level that 
varies only with time and a delta-pressure relative to a 
reference location that varies with spatial position and time; (2) 
The pressure appearing in the equation of state (for example, 
the ideal gas equation of state) is the mean pressure that varies 
only with time. Therefore, from the equation of state, density is 
a function only of the mean pressure level and the temperature 
(which varies with spatial position and time); and (3) The 
pressure appearing in the momentum equations is just the 
delta-pressure relative to the reference location. 

Thus, Fedorchenko’s paper provided the basis for 
modifying the CAST code. A technique from an early finite- 
difference one-dimensional flow non-acoustic code (Tew, 
Jefferies and Miao, 1978, Tew, 1983) was incorporated to 
ensure an accurate mass balance. It was not necessary to 
modify the incompressible form of the SIMPLE algorithm used 
in CAST since density does not vary with the delta-pressure 
from the reference location. More details are available in (Tew & 
Ibrahim, 2001) 

The modified CAST can now be used to simulate: (1) fixed- 
boundary, incompressible flow with an inlet and an outlet 
boundary, (2) incompressible flow with a moving boundary 
(piston, for example) with one inlet/outlet, or (3) compressible, 
non-acoustic, flow in a completely enclosed volume with a 
moving boundary. 

CFD-ACE+ Code 

CFD-ACE+ is a commercial CFD code product of CFD 
Research Corporation. It is a set of computer programs for 
multi-physics computational analysis. The programs provide an 
integrated geometry and grid generation module, a graphical 
user interface for preparation of the model, a computational 
solver for performing the simulation, and an interactive 
visualization program for examination and analysis of the 
simulation results. The standard CFD-ACE+ package includes 
the following applications: geometry and grid generation, 
advanced polyhedral solver and a post processor. 

CFD-ACE+ has been used for design and analysis of a 
wide variety of industrial applications. It can model 2-D or 3-D 
geometries and includes acoustics. 


ANALYSIS 

For this paper analyses were performed using CFD-ACE+ 
and the Modified CAST code for gas spring and 
piston/cyiinder/heat exchanger (or two-space) cases (See 
Appendix A for test rig dimensions originating from Kornhauser, 
1989): 

(1) Gas Spring case: The computational domain is 2-D and 
axisymmetric, with minimum dimensions of 0.0762 (x-direction) 
*0.0254(m) Piston stroke is 0.0762 m. Also 86*40 grids and 
200 time-steps/cycle were used. 

(2) Two-space case: The computational domain is 2-D and 
axisymmetric, with dimensions of -0.445 (x-direction) 
*0.0254(m). The piston stroke was also 0.0762 m for this case. 
148*48 grids and 960 time-steps/cycle were used. 

Helium was used as the working fluid, with the standard k-e, 
turbulence model and upwind finite differencing scheme. 


RESULTS 


Gas Spring Hysteresis Losses 

For a gas spring, the hysteresis loss is the work that is 
dissipated by the spring per cycle at steady operating 
conditions; it's also equal to the heat generated in and 
transferred out of the spring. A good way to compare 
computational and measured hysteresis losses is via plots of 
dimensionless work as a function of oscillating flow Peclet 
number. Dimensionless work and oscillating flow Peclet 
numbers are defined, respectively, as follows: 
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Where: 

c p specific heat at constant pressure, 

D h hydraulic diameter (=4 x wetted area^wetted perimeter), 

P mean spatial pressure, 

P a amplitude of mean spatial pressure 

p arithmetic mean of max. and min. mean spatial pressures 

V volume 

V 0 arithmetic mean of maximum and minimum volumes 

w loss non-dimensional work or hysteresis loss 

Greek letters 

y ratio of specific heats of fluid 
A molecular thermal conductivity 

p density, fluid mass per unit volume 

c o angular velocity (rad/sec) 
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Figure 1: Comparison Among CFD Codes for 
Dimensionless Losses & Kornhauser's (1989) 
Experimentally Derived Losses. 

Figure 1 shows the dimensionless work plotted versus the 
Peclet number for the computations of CFD-ACE+, Modified 
CAST and Recktenwald (1989) together with Kornhauser’s 
(1989) data. Only four data points were chosen for CFD-ACE+ 
over a range of Pe from 0.97 to 843. The CFD-ACE+ 
predictions fit very well with all the other results presented in 
Figure 1 . More detailed comparisons between the 
computational results and experimental data are given in a 
companion paper (Tew and Ibrahim 2001). 

Two-Space Test Rig Data and Calculation Comparisons 

Figures 2a and 2b show CFD results obtained from CFD- 
ACE+ and Modified CAST code respectively. The data are for 
the temperature difference between the wall and the center of 
the heat exchanger at the following different axial locations: 
entrance of the heat exchanger and 1/16, 1/8, 1/4 & 1/2 of the 
annulus total length, measured from the entrance. 



Figure2a: CFD-ACE+ Calculations of Temperature 
Difference from Heat Exchanger Center to Wall (148x48 
grids, 960 time steps/cycie). For Comparison with 
Kornhauser (1989) Exp. Data: Run #12071539, 201.7 RPM, 
1.008 MPa Mean Pressure. 



Figure 2b: Modified CAST Calculations of Temperature 

Difference from Heat Exchanger Center to Wall (34x20 
grids, 120 time steps/cycle). For Comparison with 
Kornhauser (1989) Exp. Data: Run #12071539, 201.7 RPM, 
1.008 MPa Mean Pressure. 

Upon comparing Figures 2a & 2b, the two codes CFD- 
ACE+ & Modified CAST behave in a similar way to the 
experimental data (see Tew & Ibrahim 2001) with some 
differences in the following: 1) Modified CAST predicts 
maximum temperature differences (center/wall), except at the 
entrance, larger than CFD-ACE+ and closer to the experimental 
data, 2) Also Modified CAST shows effects of the axial location 
(up to 180 Degrees-Crank-Angle) more like the experimental 
data trend than CFD-ACE+, 3) Both codes show results 
independent of the axial location beyond 230 degrees-crank- 
angle, unlike the experimental data. 

Figures 3a and 3b show CFD results obtained from CFD- 
ACE+ and the Modified CAST code, respectively. The results 
are for heat transfer per unit area at the following different 
locations in the heat exchanger: entrance of the heat 

exchanger and 1/16, 1/8, 1/4 & 1/2 of the annulus total length 
measured from the entrance. 

C«J 



z 

* 

> 

Figure 3a: CFD-ACE+ Calculations of Heat Transfer per 
Unit Area at Different Locations in the Heat Exchanger 
Relative to Entrance to Cylinder (148x48 grids, 960 time 
steps/cycle). For Comparison with Kornhauser (1989) Exp. 
Data: Run #12071539, 201.7 RPM, 1.008 MPa Mean 
Pressure. 
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Figure 3b: Modified CAST Calculations of Heat Transfer 
per Unit Area at Different Locations in the Heat Exchanger 
Relative to Entrance to Cylinder (34x20 grids, 120-time 
steps/cycle). For Comparison with Kornhauser (1989) Exp. 
Data: Run #12071539, 201.7 RPM, 1.008 MPa Mean 
Pressure. 


Again, upon comparing Figures 3a & 3b, the two codes 
CFD-ACE+ & Modified CAST behave in a similar way to the 
experimental data (see Tew & Ibrahim 2001) with some 
differences in the following: 1) Modified CAST shows the effect 
of the axial location more like the experimental data trend, 
however, 2) except at the entrance, the CFD-ACE+ heat 
transfer/unit area is closer to the experimental data than 
Modified CAST 

The results obtained at this point emphasize the fact that in 
order to obtain a reliable CFD simulation for the Stirling engine, 
there is a need for accurate experimental data. The section 
below describes the test facilities at CSU & UMN and the plan 
for running code validation experiments. 

CODE VALIDATION EXPERIMENTS 

CSU-SLRE Test Rig The CSU Stirling-Laboratory-Research- 
Engine (SLRE) derives from a two-piston Stirling-engine test rig 
developed by NASA Jet Propulsion Laboratory in 1982 (Hoehn, 
1982). A special test section is being designed to be bolted 
between the two existing piston-cylinder flanges on the rig. The 
pistons will drive atmospheric-pressure air or helium through 


the test section at Reynolds and Valensi numbers 
representative of typical Stirling engine practice. Within the test 
section will be various passages representative of regions of 
the Stirling engine where flow area changes abruptly or flow 
velocity changes direction. This may include the ports leading to 
the compression space or the transition between tubes or 
rectangular passages leading to the expansion space. In the 
test rig, the passages will be two dimensional, either axial (x, r) 
or rectangular (x, y), so that the flows can be compared directly 
with 2-D codes. 

Flow measurement will be either by LDV (Laser-Doppler 
Velocimetry) or PIV (Particle imaging Velocimetry), through 
transparent walls in the test section. Both methods require 
tracer particles within the gas flow. LDV measures individual 
tracer particle velocities one after another as it encounters them 
within its measurement volume, whereas PIV measures 
trajectories of many particles simultaneously within a 
measurement plane. PIV is the preferred method because it 
provides flow field rather than single point information. For best 
compatibility with laser measurements, test sections with flat 
walls are best. A cylindrical test section requires optical 
correction for measurements of other than the axial component 
of velocity. 

The following illustration shows a cylindrical test section 
comprising a transparent acrylic tube with an insert that 
produces an annular flow passage leading to an abrupt 
expansion at the center of the section. It will probably not be 
built, due to concerns about laser measurements through the 
circular wall. A rectangular test section has yet to be designed. 

UMN Test Rig 

The University of Minnesota test rig is an oscillatory flow 
facility, which is scaled up in size and down in frequency so the 
actual system is dynamically simulated. The drive section of 
the facility is a simple piston cylinder arrangement. The piston 
is driven by a scotch yoke arrangement so that its motion is a 
sinusoidally oscillating flow with a time mean bulk velocity of 
zero. The dimensionless amplitude and frequency match those 
of the engine through the Reynolds number based on peak 
velocity of the cycle and the Valensi number. This facility has 
been used for Stirling engine development under the SPRE 
Program (Qiu and Simon 1994). The facility is larger than the 
CSU-SLRE engine, which allows detailed measurements and 
clear visualization. It cannot match the engine frequency and 
must rely on dynamic similarity. The test section is a 
representation of a section of an engine. It consists of a 
cylindrical extension of the piston-cylinder region of the drive 
section and a radial channel between two disks, one of which is 
attached to that cylinder. It represents flow passages that could 
be found in either of the two ends of the engine, cold or hot 
section, in which the flow is driven 
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Figure 4 A Cylindrical Test Section Comprising a 
Transparent Acrylic Tube with an Insert, designed for the 
CSU-SLRE Test Rig. 

axially then radially, as it is pumped from the compression or 
expansion space radially outward toward the regenerator, then 
reverses. Flow through geometries of this type seems to not be 
computed well and there seems to be a need to improve 
turbulence closure and flow transition modeling. Measurements 
will include unsteady velocities, resolved in space and time. 
Further documentation includes visualization and measurement 
of unsteady flow separation and reattachment through the 
abrupt flow geometry change regions. Instruments which may 
be brought to bear on the problem are hot-wire anemometry, 
laser-Doppier anemometry and flow visualization with smoke, 
tufts and neutrally buoyant helium bubbles. Time records of the 
measurements will be made then decomposed according to 
time within the cycle and ensemble-averaged. 

CONCLUDING REMARKS 

Comparisons were made between the commercial CFD- 
ACE+ code and compressible non-acoustic calculations using 
Modified CAST, developed at CSU. The comparison was made 
for gas spring as well as two-space rig experimental data. The 
two codes predicted similar trends to the experimental data for, 
Work-Loss (gas spring), temperature differences (wall-center in 
the heat exchanger annulus) and wall heat flux per unit area. 
However, the two codes showed significant differences under 
certain conditions for the two-space test rig data. The reasons 
will be explored further. This work emphasizes the need for 
conducting accurate code validation experiments for oscillatory 
flow as it occurs in current Stirling engines. 
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APPENDIX 

Gas Spring & Two-Space Test Rig Dimensions 


Table Al: Gas Spring Dimensions 


Physical Quantity 

Symbol 

Value 

Cylinder Bore (Diam.) 
Piston Stroke 
Volume Ratio 

D 

S 

'V 

50.80 mm (2 in.) 
76.2 mm (3 in.) 
2.0 

Table A2: Two-Space Test Rig Dim. (Physical & Simulated) 

Physical Quantity 

Physical Value 

Simulation Value 

Cylinder Bore 
Piston Stroke 
Volume Ratio 
Annulus O.D. 
Annulus I.D. 
Annulus Gap 
Annulus Length 
Min Pist/Head Clr. 

50.80 mm (2 in) 
76.20 mm (3 in) 
2.0 

44.5 mm (1.75 in) 
39.4 mm (1.55 in) 

2.5 mm (0.10 in) 
445 mm (17.5 in) 
2.9 mm (0.11 in) 

50.80 mm (2 in) 
76.20 mm (3 in) 
2.0 

50.80 mm (2 in) 
46.4 mm (1.83 in) 
2.2 mm (0.09 in) 
445 mm (17.5 in) 
2.9 mm (0.11 in) 
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